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Abstract. 

In dusty (complex) plasmas, containing mesoscopic charged grains, the grain-grain 
interaction in many cases can be well described through a Yukawa potential. In this 
. Review we summarize the basics of the computational and theoretical approaches 

H I capable of describing many-particle Yukawa systems in the liquid and solid phases 

and discuss the properties of the dynamical density and current correlation spectra of 
^ , three- and two-dimensional strongly coupled Yukawa systems, generated by molecular 

• ^ ' dynamics simulations. We show details of the uj{k) dispersion relations for the collective 

excitations in these systems, as obtained theoretically following the quasilocalized 
[ charge approximation, as well as from the fluctuation spectra created by simulations. 



The theoretical and simulation results are also compared with those obtained in 
complex plasma experiments. 



> 

G^ ' 1. Introduction 



Strongly coupled plasmas - in which the average potential energy per particle dominates 



00 
O . 

OO , over the average kinetic energy - appear in a wide variety of physical systems: dusty 
^ ' plasmas, charged particles in cryogenic traps, condensed matter systems such as 



molten salts and liquid metals, electrons trapped on the free surface of liquid helium, 
astrophysical systems, such as the ion liquids in white dwarf interiors, neutron star 
5^ ! crusts, supernova cores, and giant planetary interiors, as well as in degenerate electron 
or hole liquids in two-dimensional or layered semiconductor nanostructures [I]. Many 
of these systems share some properties, which allows modeling them by considering 
explicitly only a single type of charged species and using a potential that accounts 
for the presence and effects of other types of species. This latter may be thought of as 
forming a charge-neutralizing background, which is either non-polarizable or polarizable. 
In the first case the interaction of the main plasma constituents can be expressed by 
the 0(r) oc 1/r Coulomb potential, while in the case of polarizable background the use 
of the 0(r) cx exp(— r/AD)/T Yukawa potential is appropriate to account for screening 
effects (Ad is the Debye length). Perhaps the most important realizations of systems 

I To whom correspondence should be addressed (donko@mail.kfki.hu) 
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lending themselves to the approximation of the interaction by the Yukawa potential are 
charged colloids [21 [3], HI |5] and dusty (complex) plasmas [6j (for comprehensive review 
on dusty plasmas see e.g. [3 E])- 

In the case of 2D colloidal systems the microscopic particles move in thin liquid 
films or between two closely separated glass plates. 

In dusty plasmas both three-dimensional (3D) and two-dimensional (2D) settings 
appear in nature and in laboratory environments. In laboratory experiments 2D systems 
appear as a layer of dust particles levitated in gaseous discharges. While most of the 
studies on this latter system have been carried out in the crystalline state (for early 
references to "plasma crystals" see P UHl [HI 112]), the liquid state is receiving more 
current attention [l3l [HI [151 [T6|, [T71 [H] . The important difference between colloidal and 
dusty plasma systems is in the damping rate in particle dynamics and concomitantly 
in the wave dynamics. In colloidal suspensions the background liquid exerts a large 
friction on the moving charged particles, while in dusty plasmas the background is 
gaseous and therefore the friction is lower and the damping of the waves is weak. For 
this reason our focus in this Review will be directed at dusty plasmas. The Review will 
cover studies of strongly coupled plasmas mostly in the liquid state, where both free 
motion and localization intervene. The principal observation is that from the point of 
view of collective behavior it is the localization - even though imperfect localization, or 
quasilocalization - of particles that plays the principal role. In contrast to the Vlasov 
plasma where the collective modes arise from the fluid-like continuum behavior, in the 
strongly coupled liquid they are more related to the normal modes of the interacting 
quasilocalized particles. This, of course, suggests a link with the harmonic phonon 
theory of crystal lattices. At the same time, one has to allow for the randomness of the 
distribution of the particles and for the finite lifetime of the localization in the constantly 
changing potential landscape. This latter process is expected to be primarily responsible 
for the damping of the collective modes, in contrast both to Vlasov plasmas, where 
Landau damping dominates and to weakly correlated plasmas where coUisional damping 
is the principal damping mechanism. This physical picture suggests a microscopic 
equation-of-motion model where the particles are trapped in local potential fluctuations. 
The particles occupy randomly located sites and undergo oscillations around them. At 
the same time, however, the site positions also change and a continuous rearrangement of 
the underlying quasi-equilibrium configuration takes place. Inherent in this description 
is the assumption that the two time scales are well separated and that for the description 
of the fast oscillating motion, the time average (converted into ensemble average) of the 
drifting quasi-equilibrium configuration is sufficient. Here the distinction between the 
"direct" and "indirect" thermal effects should be emphasized: the former are responsible 
for the actual motion and migration of the particles, the latter refer to the accessibility 
of the possible configurations of the random sites and to the temperature dependence 
of the probability of a particular configuration. 

The interaction potential energy of particles in Yukawa liquids is given by (e.g. 
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exp(-r/AD) 



Aneo r ' ^^"^ 
where Q is the charge of the particles, is the permittivity of free space, and the 
Debye length Ad accounts for the screening of the interaction by other plasma species. 
The main (dimensionless) parameters, which fully characterize the systems are: (i) the 
coupling parameter (defined in the same way as for Coulomb systems): 

1 

r = -1 — (2) 

AneoakBT' ^ ' 

where a is the Wigner-Seitz (WS) radius and T is the temperature, and (ii) the screening 
parameter 

- = f. (3) 

r is the customary measure of the ratio of the average potential energy to the average 
kinetic energy per particle; the strong coupling regime, relevant here, corresponds to 
r ^ 1. In the K — limit the interaction reduces to the Coulomb type, while at 
/€ — > CX3 it approximates the properties of a hard sphere interaction. 

At K = 0, in 3D, the liquid domain is limited to coupling parameter values F < 175 
pm [211 [22] . where the plasma is known to crystallize into a bcc lattice [231 121! • ■'■^ 2D, 
crystallization into hexagonal lattice occurs at a lower value of coupling, at F ~ 137, 
as found by computer simulations [251 ES] and proven by experiments [27] as well. At 
K > 0, 3D systems may crystallize either in a bcc or in a fee lattice, depending on k, 
as found by Hamaguchi et al. p8j in their calculation of the phase diagram of Yukawa 
systems. In 2D the crystallized form of the systems is always hexagonal. 

The possibility of characterizing a Yukawa liquid with an effective coupling 
coefficient F* instead of a (F, k) parameter pair has recently been addressed. In 3D 
F* was derived by Vaulina, Fortov and co-workers [291 EOl [31] on the basis of the 
frequency of dust lattice waves. Subsequently, for 2D Yukawa liquids the F* = /(F, h) 
relationship was established by prescribing a constant amplitude for the first peak of 
the pair correlation function for fixed values of F* [32j. The solid-liquid melting line 
(rmeiting VS. K) in both 3D and 2D system was found in these studies to follow closely 
constant F* values. 

Additional characteristic parameters of the systems investigated here are the WS 
radius a and the plasma frequency cuq, which are given for 3D and 2D systems by: 

asD = (4n3D7r/3)-i/=^ (4) 
a2D = (^2D7r)"^/^ (5) 

where nsD and n2D are the 3D (number) density and the 2D areal (number) density of 
particles, and 



t^0,3D = \\ (6j 
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and 



(7) 



'^0,2D — 



2eoma2B 



Note that in 2D the nominal plasma frequency ujo may also have different definitions, 
and some of the authors use the lattice constant instead of the WS radius as a length 
scale. 

Other important frequencies characterizing strongly coupled plasmas are the 
Einstein frequencies, which are the normal modes of oscillation of a test charge in the 
presence of a given (static) distribution of charges. Einstein frequencies are well known 
for lattice structures, however, there has been relatively little work done on disordered 
and hquid phase systems |33j. Such systems are being studied through the combination 
of theoretical and simulation approaches [SH [35] . 

As to the possibilities of theoretical description, many-body systems can be treated 
theoretically in a straightforward way in the extreme limits of both weak interaction 
and very strong interaction. In the first case, one is faced with a gaseous system, or 
a Vlasov plasma, where correlation effects can be treated perturbatively (F ^ 1). In 
the case of very strong interaction, the system crystallizes, the particles are completely 
localized and phonons are the principal excitations. In the intermediate regime - in the 
strongly coupled liquid phase - the localization of the particles in the local minima of 
the potential surface still prevails, however, due to the diffusion of the particles the time 
of localization is finite [31]. The localization of the particles (which may typically cover 
a period of several plasma oscillation cycles) serves as the basis of the Quasi- Localized 
Charge Approximation (QLCA) method [36l [37] . Besides the theoretical approaches 
computer simulations have proven to be invaluable tools for investigations of strongly 
coupled liquids of charged particles. Monte Carlo (MC) and molecular dynamics (MD) 
methods have widely been apphed in studies of the equilibrium and transport properties, 
as well as of dynamical effects and collective excitations. The main difference between 
the two techniques is that in an MC simulation the particle configuration with the 
lowest energy is searched for, whereas MD simulations provide information about the 
time-dependent phase space coordinates of the particles, this way allowing studies of 
dynamical properties. 

This paper intends to review the dynamical properties and collective behavior of 
strongly coupled Yukawa systems in the liquid and solid phases, in two and three 
dimensions. First we describe the numerical as well as theoretical methods used, in 
Sections [2] and [3l respectively. The analysis of the collective mode behavior of 3D 
liquids is presented in Section 14.11 In the 2D case we investigate both an ideally narrow 
particle layer, and a layer having a finite width, where particles are confined by an 
external parabolic potential. The analysis of these systems is described in Sections 
14.21 and 14.31 respectively. Section [5] gives a brief summary of the experimental studies 
relevant to the theoretical and simulation results reviewed. Section [6] gives a summary 
of the paper as well as a short outlook on the topics, which may have been additional 
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subjects of this Review. 

2. Molecular dynamics simulations 

Molecular dynamics simulations follow the motion of particles by integrating their 
equations of motion while accounting for the pairwise interaction of the particles, as 
well as for the forces originating from any external field(s), see e.g. |38]. In the plasma 
/ gas background environment friction forces and randomly fluctuating forces also act 
on the particles in addition to the forces arising from the interaction of the electrically 
charged particles. The general form of the equation of motion (of a "test particle i) is 



where Fj is the force originating from the interaction with particle j, Fext is the force 
originating from any external field, rj is the friction coefficient, and R represents a 
Brownian randomly fiuctuating force (Langevin force). The results presented here 
correspond to "idealized" Yukawa liquids for which 77 = and R = are assumed. 
Also, in most of our studies we investigate infinite (unconfined) systems, for which 
Fext = 0, as well, although as an example a quasi-two-dimensional liquid - confined by 
an external parabolic field - is also studied. In the case of unconfined systems periodic 
boundary conditions (PBC) are imposed in the simulations. In the case of confinement 
along one of the coordinates, PBC-s are used in the unconfined directions. It is noted 
that for charged colloids Brownian molecular dynamics simulation [3], HQ] is widely used, 
which represents the extreme limit of large friction and large R. In this case the inertial 
term (mrj) in eq. (E]) is neglected. 

The calculation of the force acting on a particle of the system, Fj, is relatively 
simple in the case of short-range potentials (like the Lenard- Jones potential or the 
potential). In this case MD methods make use of the truncation of the interaction 
potential thereby limiting the need for the summation of pairwise interactions around a 
test particle to a region of finite size. In the case of long-range interactions (e.g. Coulomb 
or low-fi: Yukawa potentials), which are also of interest here, however, such truncation of 
the potential is not allowed, and thus special techniques, like Ewald summation [H] , have 
to be used in MD simulations. Besides the Ewald summation technique there exist few 
additional methods, like the fast multipole method and the particle-particle particle- 
mesh method (PPPM, or P3M), which can be used to handle long-range interaction 
potentials, see e.g. [?2]. The results presented here for Coulomb systems have been 
obtained from simulations using this latter method [13[ HH HHl Il6[ HT] . In the PPPM 
scheme the interparticle force is partitioned into (i) a force component FpM that can be 
calculated on a mesh (the "mesh force") and (ii) a short-range ("correction") force Fpp, 
which is to be applied to closely separated pairs of particles only. In the mesh part of the 
calculation charged clouds are used instead of point-hke particles and their interaction 
is calculated on a computational mesh, taking also into account periodic images (for 



(see e.g. [39]): 




(8) 
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more details see [IHl HZ])- This way the PPPM method makes it possible to take into 
account periodic images of the system (in the PM part), without truncating the long 
range Coulomb or 1ow-k Yukawa potentials. For screening values k > 1 the PP part 
alone provides sufficient accuracy. In these cases the mesh part of the calculation is not 
used, the interaction forces are summed for particles situated within a (/t-dependent) 
cutoff radius around the test particle. Identification of these "neighboring" particles is 
aided by the "chaining mesh technique" . 

In the simulations presented here usually a spatially random particle configuration 
is set up at the initialization, with particle velocities sampled from a Maxwellian 
distribution of temperature Tq, which corresponds to the desired value of the coupling 
parameter F [see eq. ([2])]. The equations of motion of the particles are integrated using 
the leapfrog scheme or the velo city- Ver let scheme. The desired system temperature 
is reached by rescaling the particle momenta during an initialization phase of the 
simulation. In equilibrium MD simulations measurements on the system are taken 
following this phase, in the state of thermodynamic equilibrium. During this phase 
thermostation is usually no longer applied. If thermostation is necessary, rescaling of 
particle velocities is to be avoided, algorithms like the Nose-Hoover thermostat can be 
applied (see e.g. [Ml SHI iS]). 

In our studies measurements on the system are taken at constant volume {V), 
particle number (A^) and total energy (E). The MD simulations directly provide the 
pair correlation function (PCF) of the system, which is the basis for the calculation 
of thermodynamic quantities (not detailed here, see e.g. [32]), and is also required as 
input to the QLCA equations for the calculation of the dispersion relations and other 
quantities (see later). 

In the MD simulation information about the (thermally excited) collective modes 
and their dispersion is obtained from the Fourier analysis of the correlation spectra of 
the density fluctuations 

p{k,t) = k exp[ikxj{t)] (9) 
j 

yielding the dynamical structure function as [50] : 

S(k,uj) = ^^ lim -r^|p(/t,u;)P, (10) 

where AT is the length of data recording period and p{k, u) = J^[p{k, t)] is the Fourier 
transform of (Q. 

Similarly, the spectra of the longitudinal and transverse current fluctuations, 
L{k,uj) and T{k,uj), respectively, can be obtained from Fourier analysis of the 
microscopic quantities 

X{k,t) = k (t) exp [ikxj (t)] , 

j 

T{k,t) = k 'y^^Vjyjt) exp[ikxj(t)], (11) 
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where Xj and Vj are the position and velocity of the j-th particle. Here we assume that 
k is directed along the x axis (the system is isotropic) and accordingly omit the vector 
notation of the wave number. The way described above for the derivation of the spectra 
provides information for a series of wave numbers, which are multiples of /Cmin = 2tt/H, 
where H is the edge length of the simulation box. The collective modes are identified as 
peaks in the fluctuation spectra. The widths of the peaks provide additional information 
about the lifetimes of the excitations: narrow peaks correspond to longer lifetimes, while 
broad features are signals for short lived excitations. 

3. Theoretical approaches 

The Molecular Dynamics calculations compute the dynamical density-density and 
current-current correlations (dynamical structure functions), from whose behavior the 
dispersion relations for the collective modes can be inferred. Following the same route 
in a theoretical analysis would be an extremely ambitious undertaking. Calculating 
the dynamical structure functions is not an easy task and not much progress has been 
achieved so far along this line. The single-particle and collective microscopic dynamics 
of a classical 3D Yukawa fluid was first analyzed by Barrat et al. [51], on the basis of 
memory function and mode-coupling theories. They have found that the longitudinal 
current fluctuations and the velocity autocorrelation function cross over continuously 
from the behavior characteristic of classical fluids with short-range interactions to the 
dynamics of a one-component plasma as the screening parameter k of the Yukawa 
potential is reduced. 

2D Yukawa systems in the liquid phase were considered by Lowen [IQ] and Murillo 
and Gericke [52j. In this latter work radial distribution functions have been computed 
with the hypernetted chain equations and were compared with those obtained from 
molecular dynamics simulations. The dynamical structure function obtained from 
the RPA approach extended by local field corrections was shown to be inadequate to 
reproduce the features of the structure function obtained from molecular dynamics. Ref. 
[lO] focused mostly on the static properties and Brownian dynamics of the system, while 
also considering some features of the dynamical fluctuations. Applying the viscoelastic 
approximation Murillo [53] analyzed some aspects of the transverse current fluctuations. 

Fortunately, for the determination of the collective mode spectrum a much more 
direct approach, via the analysis of the dielectric response (tensor) function, is available. 
Thus the primary goals of the analytical methods discussed below are the determination 
of the dielectric function and the derivation of the ensuing dispersion relation for the 
collective modes. 

The dielectric tensor in the spatially homogeneous liquid phase is diagonal in the 
coordinate system, where k is along one of the coordinate axes. Accordingly, the 
collective modes can be classified by their polarization into longitudinal and transverse 
modes. In the crystalline solid phase the rotational symmetry is broken, the structure of 
the dielectric tensor is more intricate and the longitudinal and transverse polarizations 



Dynamical correlations and collective excitations of Yukawa liquids 



8 



do not, in general, represent eigenpolarizations anymore. In this Review we are 
concerned with the collective mode structure of the liquid phase, but the understanding 
of the behavior of the collective modes in the solid phase has a bearing, as we will 
discuss, on the formation of the collective modes in the strongly couple liquid phase as 
well. 



3.1. Fluctuation-Disspation Theorem 

The link between the ^(k, a;), L(k,u), T(k, cu) spectra measured in the simulations and 
the dielectric function is provided by the Fluctuation-Dissipation Theorem 

^(k,^) = —L{k,uj) = — — ImxLlk,^) = — - — (12) 



UJ 



^(k,u;) = 7^—5 — ImxT(k,u;), 
/c^ -Kpnuj 

where jd = l/kT, x^,^(k, c<j) is the susceptibility tensor, and x^j,(k, c<j) is the proper (or 
total) susceptibility tensor. 

At the fl value where the dispersion relation is satisfied, Xl T(k, fi) vanishes. This, 
in general, happens only at a complex frequency, the imaginary part of which being 
characteristic of the damping of the mode. Since the dynamical structure functions 
are plotted and analyzed for real frequencies only, XLT(k, f2) reaches only a minimum 
at some value of the real u, which can be expected to be in the vicinity of the actual 
complex Q: this is the frequency value that can be identified at which the peak of the 
fluctuation spectrum occurs. 



3.2. Dielectric Response Function 



The tensorial dielectric response function e^j,(k, cu) can be expressed either in terms of 
the susceptibility tensor x^uiM, ^^) or the proper (or total) susceptibility tensor XMi^(k, uj) 
and the Fourier transform (p{k) of the interaction potential ([1]). This latter depends on 
the dimensionality of the system. In 3D 

1 



and in 2D 



Then 



Eq k^ + ^2 
1 



2^0 (A;2 



K 



2^1/2 • 



(13) 



(14) 



e^^(k, Lu) = 6^^ - (p{k)x^^„(k, cu). (15) 

In the coordinate system where k is along the z axis in 3D and along the y axis in 2D, 
the isotropic hquid e^i,(k,u) has the structure 

ETiKu) 

eT{k,uj) 
£L(k,c^) 



'/lU 



3D 



(16) 
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2D 



eT(k,uj) 
e^{k,uj) 

and the dispersion relations for the collective modes are given by 

£L(k,cu) = 0, (a) (17) 
e^\k,uj) = 0. (b) 

The longitudinal dielectric function has the immediate physical significance that it 
relates the externally imposed electric field to the total (external+polarization) field by 
Etotaiik, u) = Eey^ternaiik, u) / Ehik, u) . In coutrast, the transverse dielectric function has 
well-defined physical meaning only in terms of the full electrodynamics of the system 
[5^ . Here, there is a certain degree of arbitrariness in the definition of eTO<;Uj). A 
useful alternative formulation of the dispersion relations is in terms of the external 
susceptibility 

XLk,^ = — — — , 18 
XTik,uj) = XTiKuj). 

This allows expressing the condition for the collective excitation in the universal form 
XL,Uk,^) = 0. (19) 

Xfiu{k,Lj) embodies all the dynamical properties of the system, which stem partly 
from interparticle correlations, partly from the random motion of the particles. Over 
the past half century an immense effort has gone into the calculation of this quantity 
for Coulomb systems, both classical and quantum. Most of the work focused on weakly 
coupled (r <^ 1) or moderately coupled (1 < F < 10) systems. Interest in strongly 
coupled Coulomb and Yukawa systems is more recent [55l |56] . In the strongly coupled 
domain the dynamics is dominated by correlations. Here and in the sequel we will 
mostly ignore the effect of thermal motions on Xfj.u{k,uj); some comments on how to 
abandon this simplification will be made later in this Section. 

While our focus in this Review is on the strongly coupled liquid phase, it will be 
instructive and of interest to begin with an orientation based on the weakly coupled 
Random Phase Approximation (RPA) theory. The RPA or Vlasov description is based 
on the assumption that the mean field dominates the particle-particle interaction and 
correlations can be ignored. This is tantamount to taking Xfiu{k, uj) as that of the non- 
interacting gas (although perhaps not quite obviously: for a discussion see e.g. [57]): 
Xfj,u(k.,u) = xo(k,u!)6fj,u, which, with the neglect of thermal motion is 

,^ . Ti k , > 

Xok,a;= -. (20) 

This leads to the simple expressions for the elements of the dielectric tensor 

eL(k,^)=eT(k,a;) = l- , , 3D (21) 

eLi'k,u) = ETiKu) = 1 - — , 2D 
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where co'o,3d and 0^0,20 are the respective 3D plasma frequency and the 2D nominal 
plasma frequency defined in Eqs. and ([7]); A; = ka. 

The dispersion relations (and their small-A; approximations) for the 3D and 2D 
longitudinal modes follow immediately from (fT7j) : 

J 2 2 

f^o,3D(k) = ^o,3DTi— ^ ^ 3D (22) 

J2 2 
02 _ 2 ^ ^0,2D 

^'0,201*^^ — ^0,2D n1/2 

For /c — > the longitudinal mode is acoustic, i.e. c<jL(/i; ^ 0) = sk, with the 3D and 2D 
acoustic velocities ssd and S2d: 

S3D = , (23) 

K 

S2D = 



Note that if we compare 3D and 2D systems with the same average interparticle distance, 
the acoustic speed in 2D is different by a factor y^f^ than in 3D. The acoustic behavior 
in 2D is of course at complete variance with the corresponding — > of an unscreened 
Coulomb plasma, i. e. the limit k = 0, where uj oc 

It is clear that there is no mode satisfying the transverse dispersion relation (11 7b ): 
the mean field RPA model, which is devoid of correlations, cannot support a transverse 
shear wave, since shear is a fundamentally correlational phenomenon. 



3.3. Quasi Localized Charge Approximation 

While the RPA provides a description of the weakly coupled gas, the strongly coupled 
liquid state of a Coulombic or Yukawa system requires a different approach. There 
have been various attempts over the years to calculate dispersion relations and related 
quantities for such systems. Noteworthy approaches include the high frequency sum rule 
method [50], the application of the STLS (Singwi, Tosi, Land and Sjolander) technique 
originally developed for the electron gas in metals [581159]. the memory function approach 
[601 il] and the viscoelastic model [53l [621 [631 [M] • 

In the long run, from a practical perspective most of these methods have turned 
out to be problematic. The problems that occur vary: they range from weak theoretical 
foundation through being more appropriate for static than dynamical processes, to 
resulting in an unwieldy formalism. On the other hand, a method originally proposed 
by Kalman and Golden in [65] that has become known as the Quasilocalized Charge 
Approximation (QLCA) (for a review see [Ml [37]) has led to quite a successful history of 
accomplishments. The measure of success in this context is (a) the ability to calculate 
from available static data dynamical quantities that lend themselves to comparison with 
numerical or laboratory experiments; (b) solid agreement with the outcomes of MD 
simulations; and (c) a good accord with the newly available laboratory experiments 
(still in a rather limited number) on complex plasma wave propagation. The following 
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is a concise description of the QLCA method; for more details the reader is referred to 
[36]. 

The conceptual basis for the QLCA has been a model that implies the following 
assumptions about the behavior of a strongly coupled Coulomb or Yukawa liquid: (i) 
in the potential landscape within the many-body system deep potential minima form 
that are capable of trapping (caging) charged particles; (ii) a caged charge oscillates 
with a frequency that is determined both by the local potential and the interaction 
with the other (caged) particles in their instantaneously frozen positions; (iii) the 
potential landscape changes slowly to allow the charges to execute a fair number of 
oscillations; (iv) the escape from the cages of the particles is caused by the gradual 
disintegration of the caging environment; the time scale of this process is governed by 
the coupling strength F; (v) the (time and velocity dependent) correlation between a 
selected pair of particles is well approximated by the (time and velocity independent) 
equilibrium pair correlation; (vi) the frequency spectrum calculated from the averaged 
(correlated) distribution of particles represents, in a good approximation, the average 
over the distribution of frequencies originating from the actual ensemble. 

Hypotheses (i)-(iv) have undergone careful testing by a series of MD simulation 
experiments both for Coulomb and Yukawa systems, and both for 2D and 3D 
configurations [3ll [35] , which will be discussed in Section [H The validity of hypothesis 
(v) has recently been called into question in relation to multicomponent systems. The 
short time evolution of the pair correlation function in the vicinity of a particle moving 
with respect to its environment can certainly be velocity dependent and anisotropic: it 
is now believed that it is this behavior that is responsible for some discrepancies between 
MD simulation results and QLCA predictions occurring in binary Coulomb and Yukawa 
systems. It is not believed, however, that this behavior would be problematic in a single 
component system. As to item (vi), the question of the dynamical frequency distribution 
in a liquid has received very little attention, either theoretically or experimentally (the 
record is better in relation to disordered crystals, where the problem has been posed and 
approximation schemes have been proposed, although in a language where the central 
role of the dispersion relation is obscured - see, e.g. [66]). The extension of the QLCA 
in this direction, while less than pressing, would be desirable. 

The central quantity in the QLCA is the dynamical matrix, either in three 
dimensions {V = 3) or in two dimensions {V = 2): 



which is formally similar to the eponymous quantity in the harmonic theory of lattice 
phonons and is derived from the equation of motion of properly constructed collective 
coordinates. M^^{r) = d^dyCpij) is the dipole-dipole interaction potential associated 
with 0(r). D^v{k) is a functional of the equilibrium pair correlation function (PCF) 
h{r), or of its Fourier transform h{\s). 

The longitudinal and a transverse elements of the dielectric tensor are now expressed 




(24) 
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in terms of corresponding elements of D^^{]<.) : 

.,Mk,^) = l- , y (25) 

Thus the Di^(k) and -DT(k) local field functions are the respective projections of D^,^(k) 
[36] . f2o(k) is the 3D or 2D longitudinal mode frequency, found in (l22l) . One should 
keep in mind that in spite of the universality of the expression (l24l) the explicit forms 
of the 3D and 2D eL(k, a;)-s are quite different. 

We note that the input required in the calculations is the static pair correlation 
function (PCF). In earlier works PCF-s generated by the HNC (hypernetted chain 
[67] ) technique have been used as input data of the QLCA formulae to calculate the 
dispersion relations. With the advent of computer simulation techniques it turned out 
to be both more expedient and more accurate to import simulation generated PCF-s in 
the theoretical calculations. The results of the theoretical calculations presented in this 
paper use PCF-s derived from molecular dynamics (MD) computations. 

We now can examine the dispersion relations that emerge from (12^ and fl25|) in 
conjunction with (fT7k ) and (fTTb ). We will consider both the 2D and 3D cases with 
the corresponding results for the dispersion relations displayed in figures [1] and [21 
respectively. Figure [3] will compare the sound velocities and Einstein frequencies for 
these two cases. 

Turning first to the 2D case, the longitudinal dispersion relation becomes 



Qlik) = Qlik) + DM = ^o(k) - - / dVM^,(r) [e'''-^ - l] /i(r 

m J 



(26) 



0,2D 



H / [kf^Kf^ h(f)df 



(fc2 + ^2)1/2 2 

where f = r/a, and 

A2°(x, y) = ^{{l + y + y^)[l- J,(x)] + 3 {l + y + y'/S) Mx)} . (27) 

Jo and J2 are Bessel functions of the first kind. 

In the mode frequency in ( l26i) the RPA solution and the additional correlational 
part expressed in terms of the pair correlation function h{r) are clearly separated. The 
result can, however, be transformed into an alternate form, expressed entirely in terms 
of the pair distribution function g{r) = l + h{r). By introducing the extended dynamical 
matrix C^^{k): 

Qlik) = - [ dVML(r) [e'*^-^ - l] g{r) ^ (28) 
m J 

^2 /-oo 

) — / A^^ (A;f, Kf) g{f)df. 
2 Jo 



2 

2 



^0,2D 



This result shows that the RPA contribution can be interpreted in terms of the same 
physical model as the QLCA: in this unified formulation the RPA force experienced by 
the oscillating particle is due to the mean field [h{r) = 0] only. Moreover, a reflection 
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on the origin of the "1" term in the integrand identifies it as the generator of the 
Einstein frequency, the frequency of oscillation of a single particle in the frozen immobile 
environment of the other particles (see more discussion below): 

n 



E-- I d'rM^{r)g{r)=ujQ 



1 /"°° df 

,2B7^ / ^e-- [1 + Kf + i^ff] g{f). (29) 
^ Jo " 



The k —>■ behavior of the longitudinal mode is still acoustic, but the correlations 
reduce the acoustic (sound) speed below its RPA value. For k oo the mode frequency 
approaches the Einstein frequency f^E- This limiting behavior is a remarkable feature 
of strongly coupled Coulomb and Yukawa liquids [Ml EH] . 
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Figure 1. 2D Yukawa and Coulomb liquids: QLCA longitudinal and transverse 
dispersions for specified values of the effective coupling F*. 
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(a)-(d): K = K=1 K = 2 k=3 




In contrast to the weakly coupled gas described through the RPA, the strongly 
coupled liquid supports a shear maintained transverse mode. This is reflected in the 
QLCA through the transverse dispersion relation 

Q2 (k) = DT(k) (30) 



— / 0^° (/cf, Kf) h{f)df 
^ Jo 

2 / e^"" {kr,Kf)g{f)df 2D, 
Jo 
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Figure 3. 2D and 3D Yukawa and Coulomb liquids: (a,b) QLCA sound velocities and 
(c) Einstein frequencies. 



where 

e2°(x, y) = 2^{l+y + y') [1 - Mx)] - A'^'ix, y). (31) 

The last step in fl30l) follows from a simple algebraic identity and it reflects the absence 
of a mean transverse field in the medium. 

The longitudinal and transverse dispersion curves for selected k and F values are 
displayed in figure [H The ^ behavior of the transverse mode is also acoustic, 
but the correlation maintained acoustic speed is substantially below its longitudinal 
counterpart. However, the result that the transverse mode extends all the way to = 
is spurious: the liquid is unable to support a shear wave in the uniform limit. The reason 
for this flaw is well understood: it has to be sought in the neglect of the migrational- 
diffusional damping. The introduction of a phenomenological collision frequency [70] or 
of a semi-phenomenological extension of the QLCA (studied sofar only for the Coulomb 
case - see below) [M] provide an acceptable remedy. 

For k oo the transverse mode also approaches the same Einstein frequency Q-^ 
as the longitudinal mode, as dictated by the isotropy of the liquid. 

Also shown are in figure [3] the k and F dependences of the acoustic speeds and of 
the Einstein frequency. For moderate k values the sound velocities can also be obtained 
from the semi-analytic formulae |71j : 

2\8~Ydi^^ 2 dn^J F ' 
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2 \S 2 dK^ 2 dn^J r ' ^'^ ' 

where Ec = {n/2) J (j){r)[g{r) — l]dr^ is the correlation energy per particle. 

Turning now to the 3D case, the formal results of the previous derivation can, 
mutatis mutandis, be taken over, with the understanding that the explicit forms of the 
RPA frequency ^^qC^) kernels A (kf, /tf) and G (kf, Kf) are different from their 

2D counterparts. 



2 



^o(k)=c^o'3DT^^ 3D (33) 
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(34) 



(35) 



There are now 2 degenerate transverse modes. 

The expression for the Einstein frequency is also modified: 

.2 



n 
m 



d r Mi^{r)g{r) = uj, 



K 



0,3D" 



dr re 



9[r) 



3D 



(36) 



The 3D longitudinal and transverse dispersion curves for selected k and F values 
are displayed in figure El The k and F dependences of the 3D acoustic speeds and of the 
Einstein frequency are shown together with their 2D counterparts in figure [31 Finally, 
for moderate k values the sound velocities can again be obtained from the semi-analytic 
formulae [TOl [72] 



and 
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1 + nr -\ — (nr) 
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[^7(f)-l]df (37) 
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15 



re 



1 + Kr (Kr) 

2^ ' 



[g{r) - 1] dr 



(38) 



The results of a comparison between the dispersion properties of the collective 
modes in the 2D and 3D systems (assuming the same interparticle distance) may be 
summarized as follows. 

(i) While for finite n values the qualitative behaviors of the two systems are very much 
the same, there is the well known fundamental difference in the k = Coulomb 
limit between the 2D and 3D systems as to the small k dispersion of the longitudinal 
mode: uo{k — * 0) oc \/k for 2D, but in the 3D case uj{k = 0) = ct;o,3D5 the 3D plasma 
frequency. 
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(ii) Not unrelated to this difference is the behavior of the longitudinal acoustic speeds 
at finite k values: since in 2D s oc and in 3D s oc for small k, the latter 
exceeds the former by the factor a/3/2k. 

(iii) In contrast, the transverse acoustic speeds exhibit only a mild k dependence and it 
is the 2D speed that is slightly higher than its 3D counterpart. 

(iv) A similar 2D dominance prevails for the respective Einstein frequencies that govern 
the k ^ oo behavior of the modes: in 2D the Einstein frequency assumes, for any 

a somewhat higher value than in 3D. 

While (i) and (ii) are effects originating from the basic difference caused by the long 
range behavior of the Coulomb potential in a 2D vs. a 3D geometry and are already 
reflected in the RPA description, (iii) and (iv) are correlational phenomena and they 
point at the more important role the correlations play in 2D than in 3D. 

As a closing comment, it should be re-emphasized that the QLCA ignores possible 
damping mechanisms and Doppler shift, phase mixing, etc., due to the migrational- 
diffusive motion of the particles and of velocity dispersion. A method has been recently 
proposed [731 EH ESI ES] for the extension of the QLCA to take some of the neglected 
effects into account by combining the Di^{k), D^{k) as local field factors with the Vlasov 
density-density response. In this approximation 

^ . V?(fe)Xo,L,T(k,^) 

where Xo,l,t(M,^) is the longitudinal (transverse) Vlasov density response function of 
non-interacting particles. The application of this formalism to Yukawa systems has not 
been done yet, but in an early work [68] it was shown that in a 2D Coulomb system 
the combined effect of phase mixing and Landau-damping leads to the elimination of 
oscillations in the dispersion curve. The effect of Landau damping, which is not expected 
to play a major role at high F values, is probably overestimated in this work. 



3.4- Lattice phonons 

With the caveat that the present Review addresses primarily the strongly coupled liquid 
state, it will be still useful to provide an overview of the phonon dispersion in a 2D or 
3D Yukawa crystal. Such an overview will help to understand the structure of the liquid 
state in terms of a model resembling a disordered lattice and to view the collective 
modes in the liquid as being akin to the phonon excitations in the lattice. 

The phonon dispersion is traditionally calculated in terms of the lattice dynamical 
matrix defined as 

C,.{k) = --J2 ^^M^ - [e-^^-^"--^^) - 1] , (40) 

with a summation over all the lattice sites j, keeping i fixed (r^ = 0). The resemblance 
to the extended QLCA dynamical matrix is not accidental. In contrast to the 
QLCA equivalent, however, the lattice dynamical matrix reflects the symmetry of the 
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underlying lattice and not the rotational invariance of the isotropic liquid. Nevertheless, 
a dielectric tensor can be constructed along the same line in terms of the matrix D^^i^(k), 
which is defined now as C^,i,(k) with its mean field contribution removed 

D^M = C^,(k) + - /" d^r M^,(r) [e"*^'^ - l] . (41) 



This leads to a structure analogous to (I25|) : 

■ - D(k) 



(42) 



The diagonalization of e^^i, (or of C^^) is now possible in the coordinate system of the 
eigenvectors, whose orientations, in general, do not coincide either with the direction of 
k or with the crystallographic axes. 

To find the eigenmodes one can follow the traditional method (see e. g. [771 178] ) 
of solving the secular equation 

||u;2-C^,(k)|| =0, (43) 

or continue to follow the path of working with the dielectric tensor. This latter approach 
ensures that continuity with the liquid and RPA formalism is maintained. The dispersion 
relation in terms of becomes 

k^e^uK = 0, (44) 

an obvious generalization of (|T7k). In fact, except in the degenerate isotropic case, it 
also includes the transverse relation f|T7b ). 

The 2D Yukawa system crystallizes in a triangular (hexagonal) lattice. The phonon 
spectrum was first calculated by Peeters and Wu [80], followed by Wang et al. [HI]; a 
definitive calculations of the dispersion and the polarization for all propagation angles 
are given in [79l ES]. These results are shown in figure HI The mode polarizations are 
purely longitudinal or transverse for propagation along the crystallographic axes (y? = 0° 
and 30°) only, otherwise they are mixed as shown in the figure, if is the propagation 
angle measured from the axis pointing towards the nearest neighbor. The angle G 
indicated in figure H] is the polarization angle measured with respect to the propagation 
vector k. The dispersion curves are periodic in fc, but the period is simply the reciprocal 
lattice constant only along 0° and 30°; for intermediate angles it is much longer, given 
by the formula 

An 



ko = --j=\/p^+pq + q^, (45) 
where p and q are minimal integers satisfying 

The dispersion curves of the lattice and those of the strongly coupled liquid do 
not show much resemblance. Yet, if one views the liquid as an aggregate of locally 
ordered domains whose symmetry axes are randomly distributed, then the similarity to 
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Figure 4. 2D Yukawa system: lattice dispersion curves and polarizations {k = 2) 
for different angles of propagation, ip is measured from the axis pointing towards the 
nearest neighbor, is the polarization angle measured with respect to the propagation 
vector k. Partly reproduced from Ref [73], copyright (2006) by Institute of Physics 
Pubhshing. 



the liquid dispersion should be sought in a suitably averaged dispersion of the lattice. 
The strong angular dependence of the period /cq suggests that an angular average should 
generate through phase mixing a smooth dispersion. This was carried out by projecting 
out the longitudinal and transverse components of the eigenmodes and comparing their 
respective angular averages with the longitudinal and transverse liquid modes [82]. 
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Figure 5. 2D Yukawa system: angularly averaged lattice (dashed lines) and QLCA 
(solid lines) dispersions of longitudinal and transverse modes using pair-correlation 
[g{r)] data from MD simulation at F = 360 and k = 2. Reproduced from Ref. [79] . 
copyright (2006) by Institute of Physics Publishing. 

Figure [5] shows that the agreement is quite reassuring. In principle, of course, one 
has to distinguish between the spectrum of an average of configurations and the average 
of the spectra of each of the configurations: that this observation notwithstanding the 
similarity persists can be taken as an indication that the sizes of the ordered domains 
in the liquid state are sufficiently large to diminish the effect of interaction between the 
domains. 

The 3D Yukawa system crystallizes in a bcc or a fee lattice (depending on the 
value of k). A phase diagram has been given by Hamaguchi et al. [28]. Due to the 
existence of 3 rather than 2 eigenmodes and to their dependence both on the azimuthal 
and the polar angles of propagation a much more complex phonon spectrum is expected 
than in 2D. So far no systematic published study of this spectrum seems to exist; in an 
unpublished work, however, Sullivan, Kalman and Kyrkos [83] have generated a series 
of dispersion and polarization diagrams. A sample of these, for a number of ip and 
ip angles, including propagations along the principal crystallographic axes is given in 
figure [6l (G is the polarization angle measured with respect to the propagation vector 
k, (f is the polar angle in the {x, y) plane and ip is the azimuthal angle measured from 
the z axis.) Since no averaging has been performed, their comparison with the liquid 
spectrum at the present time is difficult. 




Figure 6. 3D Yukawa system at k = 1: bcc lattice dispersion curves (left column) 
and polarizations (right column). 6 is the polarization angle measured with respect to 
the propagation vector k. Continuous line: longitudinal; dashed line: first transverse; 
dotted line: second transverse polarizations. Note that the two transverse polarizations 
may be, but in general are not degenerate. <^ is the polar angle in the {x, y) plane and 
ip is the azimuthal angle measured from the z axis. 
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3. 5. Einstein frequencies 

In addition to the collective excitations, Einstein frequencies represent a dynamical 
manifestation of the strong interaction in Yukawa systems. Einstein frequencies, as 
noted above, are the frequencies of oscillation of a single particle of the system (the 
"test particle") around its equilibrium position in the immobilized frozen environment 
of the other particles of the system. For obvious reasons, from the experimental point 
of view the "freezing" of the system but one particle is not a realistic proposition. 

Thus until the advent of dusty plasma experiments Einstein frequencies were 
considered more of a theoretical construct than an observable quantity. The realization, 
however, that in the strongly coupled liquid state (but not in the crystalline sohd) they 
represent the asymptotic k —>■ oo limit of the mode dispersion has promoted the Einstein 
frequency to the rank of observable quantities [HH 185] . 

In the crystalline solid state, where the test particle occupies a lattice site, the 
assumption that the potential experienced by the test particle is a quadratic function 
of the coordinates with a positive definite second derivative is in accord with the basic 
model of the harmonic theory of phonons. The maximum number of eigenfrequencies 
of oscillation is equal to the dimensionality of the system {V = & or V = 3); because of 
the lattice symmetry induced degeneracy the actual number may be less than V. In a 
disordered lattice the degeneracy is removed and the frequencies depend on the actual 
realization of the disorder. In this case one has to distinguish between the "microscopic" 
Einstein frequencies (cue) each of which is generated by a particular realization of the 
disorder and characterized by a distribution over the ensemble, and their ensemble 
average f^E = \/ {^e)- ^^^^ latter that will be continued to be referred to as 

"Einstein frequency" in the rest of this paper. 

In addition, it is useful to consider the quantity 

V 



i.e. the sum of the squared eigenfrequencies in a particular realization. Obviously, 
(cul) = (lj^) but the distributions of and can be quite different. 

In a strongly coupled liquid the very notion of "equilibrium position" is questionable. 
Nevertheless, the "quasilocalization" condition, the basic tenet of the QLCA, is well 
satisfied for high F values, as demonstrated by MD experiments [3l] the details of which 
will be discussed below. It is in this sense that the notion of the Einstein frequency and 
its distribution can be extended to the case of the strongly coupled liquid. 

In a 3D Coulomb crystal the Einstein frequency is determined solely by the 
background, unaffected by the distribution of the (frozen) particles. This is the 
consequence of Gauss Theorem which, in turn, follows from the Poisson Equation that 
the 3D Coulomb potential satisfies. In this case 




(47) 




1 Q^n _ 1 2 



(48) 
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In a disordered lattice or in a liquid 
weaker statement 



UJt 



2^0,3D5 



is not valid anymore; it is replaced by the 



(49) 



the so-called Kohn Sum Rule [86], which also follows from the Poisson Equation. 

Thus, in a disordered system while u]^ has a spread, u]^ does not. As to the average, 
(149|) is of course also tantamount to 
1 o 



0,3D- 



(50) 



For a genuine Yukawa potential the situation is quite different. The Yukawa 
potential satisfies the screened Poisson Equation rather than the Poisson Equation. 
A useful statement can be made now only for Qe, which now can be expressed in terms 
the average of the Yukawa potential (0) as experienced by the test particle at r = 



3m 



(</>) 



(51) 
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0,3D" 



dr re 



o,3d: 
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9{r) 



dr re ^'^h{r) 



flSTj) is in agreement with fl36|) . the result obtained from the QLCA. The third line clearly 
shows that, remarkably, in the k, —>■ Coulomb limit the Yukawa Einstein frequency 
reduces to the background induced ( l50i) . even though the Yukawa system exists without 
any background. It can also be noted that = E^^t is the interaction energy density 
of the system (with [positive] Hartree plus [negative] correlation contributions). Since 
the energy is the lowest in the ordered state, the Einstein frequency must increase with 
increasing disorder. According to the known phase diagram of the 3D Yukawa system 
[28] - as already mentioned - the system crystallizes in a bcc or a fee lattice. The 
corresponding Einstein frequencies [87| 
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(52) 



constitute an absolute lower bound. 

In the 2D Coulomb system Gauss Theorem does not apply, the background plays no 
role and neither the Poisson Equation nor its screened variant is satisfied. Consequently, 
the Einstein frequency is determined by the distribution of the surrounding particles, 
both for Yukawa and Coulomb systems. In general 



"e 
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g{f) (53) 



in agreement with the QLCA result ([29 
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An argument similar to the one discussed in relation to the 3D case leads to the 
conclusion that here also the ordered state exhibits the lowest Einstein frequency. The 
lattice structure is now hexagonal, for which 

nl{K = 0) = 0.39925 ujI^^ (54) 
nl{K = 1) = 0.34433 
nl{K = 2) = 0.24347 tol^D 

These values constitute then the lowest bound for the 2D Einstein frequencies. 

In addition to the frequencies, the Einstein oscillations are also characterized by 
their eigenpolarizations. It is the distribution of the polarization angles which is of 
interest; this question has been investigated, however, only for the 2D case [82]. In 
the perfect hexagonal lattice the degeneracy of the eigenmodes renders this distribution 
isotropic. It is also isotropic in the liquid phase. However, in the intermediate range 
where the lattice disorder develops the degeneracy for the microscopic eigenmodes is 
removed and the rotational invariance of the distribution is reduced to the sixfold 
symmetry of the underlying lattice. More will be shown about this remarkable effect in 
Section 1121 

4. Simulation results 

In this Section we review the results of the extensive MD simulation work carried out 
since the beginning of this decade on the dynamical properties of Yukawa liquids. Most 
of the work was motivated by the QLCA theory and accordingly a great portion of the 
results pertaining to the areas where QLCA predictions are available are accompanied by 
comparisons with the theoretical predictions. However, the information generated by the 
simulations goes well beyond those predictions: this is eminently true for the frequency 
spectra of the dynamical density- density and current- current correlation functions 
[dynamical structure functions S{k, u), L{k, u), T{k, u)]. Beyond predictions pertaining 
to the peak positions of the spectra, identified as the frequencies of the collective 
excitations, the QLCA does not provide, apart from some qualitative estimates, any 
basis for comparison in this respect. While other works, based mostly on the memory 
function formalism [5T1 [80t ISTl [88] , have presented theoretical descriptions of some of 
the features of the structure functions, we have made no attempt to relate to these, 
rather scant, results for the purpose of comparison with simulations. 

As noted in Section 13.31 the basic hypotheses (i)-(iv) of the QLCA theory have 
undergone careful testing by a series of MD simulation experiments both for Coulomb 
and Yukawa systems and both for 2D and 3D configurations [3H [35] . With increasing 
r values, a visual inspection of the potential landscape clearly indicates the formation 
of potential wells [35]. Examination of the phase space trajectories reveals a clear 
morphological difference between low F and high F situations: in the first case the 
trajectories are open, interrupted by propagating oscillatory portions, while in the 
second case the trajectories are mostly closed and exhibit a loop structure characteristic 
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of localized oscillatory motion [31]. An example of this behavior is illustrated for a 3D 
Coulomb liquid in figures [7|(a) and (b) for F = 2.5 and F = 160, respectively. 




1 10 100 1000 1 10 100 1000 

r r 



Figure 7. Phase space trajectory segments of a test particle in a 3D Coulomb liquid at 
(a) r = 2.5 and (b) T = 160. H is the edge length of the simulation box. Decorrelation 
time of the cages (Tdocorr — i-^o^docorr/STr) as a function of T for the (c) 3D and (d) 
2D systems, for a series of k values, obtained from MD simulations. (a,b) Reproduced 
from Ref. [35 • Copyright (2002) by the American Physical Society. (c,d) Reprinted 
with permission from [Z. Donko, P. Hartmann, and G. J. Kalman, Phys. Plasmas 10, 
(5), 1563 (2003)]. Copyright (2003) by the American Institute of Physics, Ref. [35] . 

The quantification of the relationship between localization and the strength of the 
coupling has been carried out by invoking a technique due to Rabani et al. [H9] . Here a 
"cage correlation function" was introduced to characterize the gradual disintegration of 
the cage of the nearest neighbors and the escape of the caged particle. The main results 
shown in figures [7] (c) and (d) for 3D and 2D Coulomb and Yukawa systems illustrate 
the duration (in terms of plasma oscillation cycles) of the caging (decorrelation time, 
Tdecorr) as a fuuctiou of F and n. In the case of the 3D system, at k = and F = 160 
the cages decorrelate during ^ 50 plasma oscillation cycles. The decorrelation time is 
reduced to a single cycle at F 7. In the case of the 2D system it takes about 100 cycles 
for the cages to decorrelate at k = and F = 120, and we reach Tjecorr = 1 at F ^ 2.5. 
In the high-F domain we observe a strong dependence of the decorrelation time on k, 
both in 3D and 2D systems. At low values of F, however, Tdecorr depends only slightly 
on K. The decrease of the decorrelation time for increasing k can be compensated by 
increasing F, as it can be seen in figure [71^ c) and (d) [35]. It is noted that the data 
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Figure 8. (Color online) 3D Yukawa and Coulomb liquids: spectral decomposition 
of the longitudinal and transverse current fluctuations in Coulomb (a,b) and Yukawa 
(c,d) plasmas, at F = 160, k = 0, and F 380, k — 2, respectively. (The color coding 
of the amplitude is logarithmic, it only intends to illustrate qualitative features.) 



shown in figure [7] convey information about the "average behavior" of the particles, it 
is however, recognized [35] that the surrounding of individual particles may change in 
a different way, due to e.g. avalanche type excitation and migration [90]. Finally we 
note that the caging of the particles at high F values determines many of the systems 
properties as it has been discussed by Daligault for 3D Coulomb liquids [9T| . 



4.I. Three-dimensional Yukawa liquids 

The first molecular dynamics simulations on the wave dispersion relations in the fluid 
phase of 3D Yukawa systems were reported by Hamaguchi and Ohta [921 [93]. Their 
results confirmed the earlier theoretical predictions of Rosenberg and Kalman [72] on 
the longitudinal wave dispersion and were mostly in agreement with the simultaneously 
published full QLCA calculations of Kalman et al. [70]. They also demonstrated that 
the transverse wave dispersion has a cutoff at a long wavelength even in the case of weak 
screening. 

This work was followed by a series of MD simulation for the collective excitations 
in 3D Yukawa liquids to provide further comparison with the predictions of the QLCA 
theory. The simulations - of which the results are presented here for the first time - 
have been carried out using N = 12 800 - 15 625 particles. 
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To illustrate qualitatively the features of the behavior of the collectives excitations 
the spectral decomposition of the longitudinal and transverse current fluctuations is 
plotted in figure M for 3- dimensional Coulomb and Yukawa liquids. In the case of the 
Coulomb plasma, at low wave numbers the frequency of the longitudinal (£) mode is 
concentrated within a narrow frequency range [see figurelHI^a)] near the plasma frequency. 
With increasing wave number the frequency of the mode gradually spreads over a wider 
domain and shows a slightly decreasing tendency. In sharp contrast with this behavior 
the T mode frequency is spread over a wide domain, as illustrated in figure [H](b). The 
C mode of the Yukawa system is quite different from that in the Coulomb case, the 
wave frequency approaches zero at — ^ wave number. The frequency increases with 
increasing wave number up to about k = 2.0, and then starts to decrease slightly. 
Meanwhile the frequency distribution gets gradually wider. The T mode in the Yukawa 
case appears to be similar to the corresponding mode in the Coulomb system, although 
the frequency is lower, due to the weaker interaction of the particles, as a consequence 
of the screened potential. 

For a better quantitative analysis representative dynamical structure functions 
(density fluctuation spectra) S{k, uj) and spectra of the longitudinal and transverse 
current fluctuations, L{k,u!) and T{k,u), are plotted in figures [9] and [TOl respectively, 
for a high-F and a medium-F case. The S{k, u) obtained for the Coulomb case (F = 
160, K, = 0, see figure [9]^a)) peaks at nearly the same frequency for the different values of 
the wave numbers plotted, which are multiples of kmin = 0.167 (determined by the size 
of the simulation box). In the presence of screening (Yukawa potential), as shown in 
figure [9](d), the behavior of S{k,uj) changes significantly: at A; ^ the wave frequency 
tu/tuo — > [luq is defined by (E])]. The contrast between the k = and the k > cases is 
also well seen in figure [TlTa). where the dispersion curves derived from the fluctuation 
spectra are displayed. The dispersion curves for /? > are quasi-acoustic (co/ujq oc fc^^^), 
with a linear portion near k = 0, which gradually extends when n is increased. The 
(F,k) pairs for which the dispersion graphs are plotted in figure [11] have been selected 
to represent a constant effective coupling F* = 160. This definition of F* relies on the 
constancy of the first peak amplitude of the pair correlation function g{f), similarly to 
the case of 2D Yukawa liquids p32] . 




Peaks in the spectra of the compressional C mode [plotted in panels (b) and (e) of 
figures M and [10] appear at the same frequency as those in the corresponding S{k, ur) 
functions, as these functions are linked via the relation 



Compared to those characterizing the C mode, peaks in the T mode spectra are 
rather broad, as it can be seen in panels (c) and (f) of figures [^ and [TUl In the case of this 
mode there is no dramatic change between the behavior when k changes from zero to a 
nonzero value, only the mode frequency decreases, as can be observed in figure [TTT b). 

Comparison of the dispersion relations obtained from the MD simulations [via 
S{k,u)] and QLCA calculations [see equations (!33l) -( !35l) ] is presented in figure [TTl Here, 




L{k, u) 



S{k, u). 



(55) 
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Figure 9. 3D Yukawa and Coulomb liquids: density [S{k,uj)] and current [L{k,uj) and 
T{k, co)] fluctuation spectra of Coulomb F = 160, «; = (a,b,c) and Yukawa F = 200, 
K = 1 (d.c.f) systems. The curves are plotted for multiples of the smallest accessible 
wave number fcmin = 0.167. (g) and (h) show the dependence of L[k,Lu) and T[k,u)), 
respectively, on k at fixed wave number k = 1.00. (F = 360 for k = 2, and F = 1050 
for K = 3). 
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Figure 10. 3D Yukawa and Coulomb liquids: density [S{k,Lo)] and current [L{k,LLj) 
and T{k, uj)] fluctuation spectra of Coulomb F = 20, k = (a,b,c) and Yukawa F = 48, 
K = 2 (d,e,f) systems. The curves are plotted for multiples of the smallest accessible 
wave number fci„i„ = 0.156. (g) and (h) show the dependence of L{k,ui) and T{k,oj), 
respectively, on k, at fixed wave number k = 1.56. (F = 25 for k = 1, and F = 114 for 
K = 3). 
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Figure 11. 3D Yukawa and Coulomb liquids: dispersion relations for the (a) 
longitudinal and (b) transverse modes. O : r = 160, k = (Coulomb case), • 
: r = 200, K = 1, ■ : r = 380, k = 2, and A : T ^ 910, k = 3. Symbols represent 
molecular dynamics results, while the lines correspond to the predictions of the QLCA 
theory. 




K K 

Figure 12. 3D Yukawa and Coulomb liquids: (a) sound velocities and (b) Einstein 
frequency as derived from the QLCA theory (circles) and Einstein frequencies of fee 
lattice (stars) [94] . 

in the calculations of the QLCA results, we have made use of the g{r) functions obtained 
from the MD simulation. The agreement between the two sets of data is excellent for 
the C mode, while some difference in the frequency of the T waves can be seen in 
figure [TlTb). This latter may originate from the inaccurate determination of the peak 
positions of the rather broad T{k,u) spectra. It should be noted though that while the 
theoretical calculations provide an oscillatory dispersion curve for > 3 (see figure [2]), 
simulations provide reliable results (for collective excitations) for typical liquid-phase 
conditions for < 3. (At higher k values the thermal contribution in S{k,uj) apparently 
masks the collective mode peak.). The simulation results here resemble the measured 
2D dispersion curves in the liquid phase [15]. Another difference is the cutoff of the T 
mode dispersion curve at finite wave numbers. This disappearance of the shear modes 
for — > is a well known feature of the liquid state [Sni [951 EB], and the sharp cut-off 
^ for a finite k has also been observed in simulations of Yukawa systems [531 ES] ■ 
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It has been already noted that this cutoff is not accounted for by the QLCA, as it does 
not include damping effects. 

The sound velocities, derived in ( 1371) . are plotted in figure [T2l(a). while figure fT2](b) 
displays the Einstein frequency, which is defined in (l36l) . In the /t — limit (!36ll 
gives lue = ujq/\/3, and ue decreases with increasing k. For comparison, the Einstein 
frequency data of Ohta and Hamaguchi for a fee lattice [91] are also plotted in 
figure [TWb). We find an excellent agreement between the two sets of data. 
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Figure 13. 3D Yukawa and Coulomb liquids: Einstein frequency distributions for 
r* = 120: (a) r = 120, k = 0; (h) T = 150, k = 1; (c) T = 300, k = 2; (d) T = 725, 
K = 3; the vertical bars indicate values obtained by the QLCA theory, (a) also shows 
the distribution of frequencies at the lower coupling value F = 20. (e) Distribution of 
Lj^ for the same systems. 



Numerical experiments were also performed to determine the distribution of the 
microscopic Einstein frequency uje- To accomplish this, frequency histograms based on 
a few hundred, temporally uncorrelated particle configurations have been constructed. 
For the raw (particle position) data the harmonic matrix for every particle has to be 
generated: 

3+1. 
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where r^*^ is the equihbrium position of the z-th particle (local minimum of the potential 
surface), (j){r) is the interaction potential, a and /5 represent the Cartesian coordinates. 
The eigenvalues of Hais/m are the squared Einstein frequencies (3 for every particle), 
while the eigenvectors provide the polarization of the oscillation. 

A series of frequency histograms for an effective coupling parameter T* = 120 and 
different values of k are shown in figures [T3r a)-(d). The frequency distributions exhibit 
three peaks (although this is less visible in the k = 3 case). With increasing screening 
the distribution of frequencies becomes wider and its mean value is shifted towards 
lower frequency. The QLCA results for the Einstein frequency [obtained from Eq. (!36!l 
using pair correlation functions generated in the MD simulation], corresponding to the 
different values of k are also indicated in figures [T3l(a)-(d). The values are in good 
agreement with the simulation results. The effect of F at fixed (k = 0) screening is 
illustrated in figure [T3]^a). A six time decrease of the coupling parameter results in 
approximately doubled width of the Einstein frequency distribution. 

Figure [T3r e) shows the histograms for ui^, sums of the 3 microscopic squared 
Einstein frequencies, for different values of k: there is a qualitative difference between 
the K = Coulomb case where there is only a single frequency (a narrow peak) and 
the K > cases, where a distribution of frequencies is apparent. The reason for this 
difference has been discussed in Section |3l 

Further information on the collective behavior is contained in the velocity 
autocorrelation function (VACF) 

^ ^ " (|v(0)P) ' ^^^^ 
where the average is taken over the particles and different initial times. 

The behavior of the velocity autocorrelation functions of 3D Yukawa liquids 
obtained at several values of the F and k parameters is illustrated in figure [T4l Analyzing 
the behavior of Z{t) at constant n, we find a transition from monotonically decreasing 
Z{t) into an oscillating type when F is increased [figure fT4r a)]. Similarly, the shape of 
Z{t) changes drastically when k is varied at constant F, as shown in figure [T^ b). For 
more detailed analysis see [M] . 

Using the Einstein frequency u-e for the normalization of time, instead of the plasma 
frequency uq (as in figure [T4l) the Z{t) functions belonging to the same F* = 160 for a 
series of k values are displayed in figure fT5l(a). Using this normalization of the timescale 
the Z{t) curves exhibit a nearly universal behavior, where at least the first few peaks 
of Z{t) nearly overlap. This observation emphasizes the importance of the Einstein 
frequency in the dynamical behavior of the system. 

The marked oscillations of the Z{t) function in the Coulomb case is indeed expected 
on the basis of the possible coupling between the single particle motion and long- 
wavelength plasmons, whose frequencies are almost independent of k. For k > 0, 
however, even though uj{k 0) oc k, the oscillations persist; this can be explained by 
the fact that the uj{k) dispersion curve flattens at higher wave numbers, a corresponding 
peak in the frequency distribution develops. 
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Figure 14. 3D Yukawa and Coulomb liquids: (a) velocity autocorrelation functions 
at K = 1.0 and a series of T values; (b) at constant F = 100, for a series of k values. 
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Figure 15. 3D Yukawa and Coulomb liquids: (a) velocity autocorrelation functions 
for a series of k values; the time is normalized by the Einstein frequency lue- (b) 
corresponding Fourier transforms Z(lo\ for (F, k) pairs as indicated in (a). 



The Fourier transforms ^(cu) of the VACF functions (obtained at different k values 
but for constant effective couphng F* = 160) are portrayed in figure [TST b). The dominant 
peaks in the spectra - shifting towards lower frequencies with increasing k - correspond 
to the high frequency oscillations of the functions (easily observed visually). As 
discussed in previous studies (see e.g. [OH [96]) these peaks are related to longitudinal 
current fluctuations, while the broad features at low frequencies are connected to the 
transverse current fluctuations and are related to diffusion properties of the system. 
Taking the case of k = 2 as an example, figure [8](c) indicates that most of the energy of 
the L mode is concentrated around frequencies LojiOQ = 0.5, in correspondence with the 
peak of the Ziuj) function shown in figure [T5l(b). The T(fc, uj) spectra [see figure [Ht^c)] for 
any k are broader, compared to the L{k,u) spectra: the fluctuations in the transverse 
currents are distributed over a rather broad frequency domain, again in agreement with 
the behavior of the corresponding Z{uj) function. The observed features of Z{u) indicate 
an appreciable coupling between single particle motion and collective excitations in the 
2D system. 

In addition to the peaks in the frequency spectrum of the dynamical structure 
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functions a wealth of further physical information is contained in the detailed structures 
of these quantities. Of great importance would be the understanding of the evolution 
of the width of the frequency spectra as functions of k and F, since it is related to the 
damping of the collective modes. To illustrate the behavior of the dynamical structure 
function figure [16] shows the widths of the collective mode peaks as a function of wave 
number, for effective coupling values F* = 20 and 120, for k, = and 2. It is noted that 
Murillo has provided a formula for the width of the transverse current spectrum |53j . 
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4-2. Two-Dimensional Yukawa liquids 

Most of the available experimental evidence on waves in complex plasmas relates to 2D 
systems (see section [5]); much less information can be culled from observations on 3D 
systems. Thus the understanding of the collective mode structure in the different phases 
of the 2D Yukawa system has been of great current interest: over the past few years 
a substantial amount of simulation work has been performed on 2D Yukawa liquids 
[7T| [97] . In the following we present these MD simulation results on the dispersion 
properties of the liquid state and compare them with the theoretical predictions of the 
QLCA analysis of the collective modes. 

Representative density fluctuation spectra, as well as longitudinal and transverse 
current fluctuation spectra of the 2D Yukawa liquid are displayed in figure [T71 The 
dispersion curves derived from the simulation spectra S{k, u) for both modes are 
displayed in figure [181 The results shown in the latter figure at k = reproduce 
the known 2D Coulomb dispersion [95l [98]. With increasing k the mode frequencies 
rapidly diminish and the dispersion deviates more substantially from its RPA value. 
In the k —>■ limit both modes exhibit an acoustic behavior, with longitudinal and 
transverse sound velocities sl and st [SHI EH], see Eq. ( l32l) . For the longitudinal 
mode the simulation data well corroborate the theoretical predictions with the proviso 
already noted in relation to the 3D liquid: that while the theoretical calculations provide 
an oscillatory dispersion curve for k > 3 (see figure [1]), simulations provide reliable 
results (for collective excitations) for the given conditions for A; ^ 3. In the case of 
the transverse mode, the agreement between theory and MD data for moderately high 
k values is fairly good; for A; the agreement is marred by the QLCA's inability 
to account for diffusional and other damping effects [9H] that preclude the existence 
of long wavelength shear waves in the liquid state. As a result of this damping, a 
cutoff at a finite kc and zero frequency develops (a similar phenomenon was observed 
in the 3D case [531 [93]). The kc value is related to the diffusional- migrational time [98] 
through tdm = ^/K^t, where st is the transverse sound velocity. Incorporating tdm, 
calculated with the aid of the theoretically predicted st values, in the QLCA equations 
as a phenomenological damping u = I/tdm (by the u ^ u + iu replacement), good 
agreement between the theory and the MD data was restored [Tlj. The simulations 
show that the longitudinal mode is not affected by this damping mechanism: this may 
indicate that its characteristic damping time is substantially longer. 

The sound velocities and the Einstein frequency - given as the k —>■ oo limit of 
Eqs. f[26l) or f[28|) - are shown in figure [191 For comparison, also displayed is the 
thermodynamic sound velocity generated from the equation of state of a Yukawa liquid 
[32] . The sound velocities obtained here are extremely close to those of the hexagonal 
crystal [HD]- The Einstein frequency diminishes rapidly with increasing k, similarly to 
the 3D case [70l[9l]. 

It is of interest to follow the evolution of the mode structure across the liquid- 
solid phase boundary, as the isotropic liquid dispersion transits into the anisotropic 
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Figure 17. 2D Yukawa and Coulomb liquids: density [S{k,u!)] and current [L{k,Lu) 
and T{k,ui)] fluctuation spectra of Coulomb F = 120, k = (a,b,c) and Yukawa F 
= 160, K = 1 (d,e,f) systems. The curves are plotted for multiples of the smallest 
accessible wave number kmin = 0.0886. (g) and (h) show the dependence of L{k,uj) 
and T{k,uo), respectively, on k, at fixed wave number k = 1.063. (F = 360 for k = 2, 
and F = 1050 for k = 3). 



Dynamical correlations and collective excitations of Yukawa liquids 



37 




Figure 18. 2D Yukawa and Coulomb liquids: dispersion curves for (a) the longitudinal 
(£) and (b) transverse (T) modes at F* = 120 and k = 0, 1, 2, 3. Continuous 
curves: QLCA calculations; symbols: MD simulation; dashed lines: RPA dispersions. 
Reproduced from Ref. [7T]. Copyright (2004) by the American Physical Society. 



3 
CD 



1.0 
0.8 
0.6 
0.4 I- 
0.2 
0.0 



\ ' \ \ 


. ' 1 ' 1 

\ . 

\ ■ 
\ 

4 \ 

K \ 

\\ ^ 

\ 


T — 1 — 1 — r 

(a) ■ 


: 










0.0 0.5 1.0 1.5 2.0 2.5 3.0 
K 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 
K 



Figure 19. 2D Yukawa and Coulomb liquids: (a) sound velocities (heavy lines with 
circles: calculated from QLCA at F* = 120, thin lines: hexagonal crystal lattice [80] . 
dotted line: RPA values for sl, and dashed lines: 3D values at F = 160, open triangles: 
thermodynamic sound velocity) [TQl [71] ; (b) calculated Einstein frequency as obtained 
from the QLCA formula for F* = 120. Reproduced from Ref. [71]. Copyright (2004) 
by the American Physical Society. 

dispersion of the solid state. This is illustrated for the k = 2 case in figure [20l T = 500 
represents a relatively high temperature solid, where lattice defects may already show 
up, but the overall behavior (sharp separation of the mode frequencies along the x and 
y directions; compare e.g. the curves labeled Tx and Ty) reflects the conservation of the 
triangular crystalline structure. The F = 405 case corresponds to a temperature slightly 
higher than the melting temperature (our results indicate that the transition occurs at 
r = 415 for K = 2 [82]), where all long range order in the system has already been 
extinguished, but locally most of the particles sit in the somewhat distorted hexagonal 
environment. The "oscillatory" feature in the T mode around k = 2.5 can be taken as 
an indication for the transition from the ordered lattice to the disordered liquid state 
through the formation of disoriented domains of local hexagonal order. The orientation 
of these domains becomes more uncorrelated with increasing temperature. The F = 200 
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Figure 20. 2D Yukawa systems: Comparison of MD (L and T) dispersions in the 
solid phase (F — 500), just below the melting transition (F — 405) and in the liquid 
phase (F = 200) for k = 2. Shown are both x and y polarizations, where x is in the 
direction to the nearest neighbor in the hexagonal lattice. Reproduced from Ref. [52] ■ 
Copyright (2007) by the Institute of Electrical and Electronics Engineers. 



system is a typical strongly coupled liquid. Most prominent features are the isotropy 
of the dispersion (x and y directions are equivalent), and the appearance of a finite 
wavenumber cut-off for the T mode. 

The behavior of the velocity autocorrelation function Z{t) of the 2D liquid is very 
similar to its 3D counterpart, at least for short times. Representative Z{t) functions 
obtained at different screening parameter values are shown in figure [2]T a). These 
functions, when plotted against exhibit nearly universal behavior [see figure [2T](b)]. 
indicating the relevance of the Einstein frequency in determining the single particle 
properties [32]. While the in-depth analysis of the long-time behavior of the velocity 
autocorrelation function is beyond the scope of this paper, it is noted that in low- 
dimensional systems Z{t) may exhibit a slow power law decay, which makes it non- 
integrable [99]. As a consequence the diffusion coefficient may not exist for some 2D 
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Figure 21. 2D Yukawa and Coulomb liquids: (a) velocity autocorrelation functions 
for a series of k values, (b) The same data as a function of we^- The (F, k) pairs 
correspond to the same effective coupling F*. Reproduced from Ref. [32] ■ Copyright 
(2005) by the American Physical Society. 




Figure 22. 2D Yukawa and Coulomb liquids: Einstein frequency distributions for 
F* = 120: (a) F = 120, k = 0; (b) F = 160, k = 1; (c) F = 360, k = 2; (d) F = 1050, 
K — 3; the vertical bars indicate values obtained by the QLCA theory, (a) also shows 
the frequency distribution obtained at a lower coupling value F = 10. (e) Distribution 
of for the same systems. 
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Figure 23. 2D Yukawa systems: distribution of the polarization angle for the higher 
frequency normal mode for different values of the coupling parameter F across the 
crystallization boundary, k = 2. Reproduced from Ref. [104] . Copyright (2007) by 
the Institute of Electrical and Electronics Engineers. 

systems. The case of 2D Yukawa liquids has attracted considerable attention during the 
last years |100[ llOll 11021 1103] . These studies have found very nearly Z{t) oc t~^ decay 
of the velocity autocorrelation function and superdiffusion to exist for some conditions. 

Similarly to the 3D case, numerical experiments were also performed for the 2D 
case to determine the distribution of the microscopic Einstein frequencies. A series of 
frequency histograms for T* = 120 at different values of k are shown in figures [221(a)- (d) . 
We observe two peaks in the distributions, which gradually get wider with increasing 
K. The QLCA results are again in very good agreement with the mean values of 
the distributions. A wider frequency distribution appears when F is lowered [see 
figure [22r e)]. The distributions of the sums of the 2 microscopic squared Einstein 
frequencies uj^ - as shown in figure [22]^e) - are in contrast with the 3D situation. Here 
we do not find qualitative difference between the k, = Coulomb and the n Yukawa 
cases for reasons discussed in section El 

The angular distribution of the polarization vector of the higher frequency normal 
mode oscillation has also been analyzed, as an indicator of the prevailing disorder [104j . 
Figure [23] shows the distribution of the polarization angle for the higher frequency 
normal mode for different values across the crystallization boundary (at F = 415). As 
discussed in section [3l both the liquid (away from the phase transition boundary) and 
the perfect lattice (extremely high F values) exhibit a full rotational symmetry, while 
in between the sixfold symmetry of the lattice prevails. 

The diagram showing the widths of the spectra of the dynamical structure functions 
is displayed in figure [2H Comparison with data for the 3D Yukawa liquid reveals that 
the trends and orders-of- magnitudes in the two cases are not substantially different. The 
dependence of the width of the peaks in the S{k, u) spectra as a function of the reduced 
coupling parameter follows the form Alu/luq = 0.76(F*)~^/^, as it is shown in figure [251 
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Figure 24. 2D Yukawa and Coulomb liquids: spectral line width for k = and 
K = 2 at r* = 40 (a,c,e) and T* = 120 (b,d,f). Shown are FWHM (full-width at half- 
maximum) values for the most prominent peaks in the S{k,uj), L(k,uj) and T(k,u!) 
spectra. 



This monotonically decreasing function of the couphng may change character at lower 
couphng values, for which data are at present not available. This is expected on the basis 
of the prediction by Hansen et al. |50j (for the 3D case): for weak coupling the width 
is expected to increase from its Vlasov value where S{k, uo) should be extremely sharp, 
since higher coupling leads to higher collision frequency and thus to stronger damping. 
Once, however, localization sets on, further increase in the coupling is expected to 
create better localization and thus a reduction in the collision frequency and in the 
width. Thus, generation of data for lower F values would be desirable, to see whether a 
turnaround point really exists. 
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Figure 25. 2D Yukawa and Coulomb liquids: line widths of the S{k,uj) spectra as a 
function of the reduced coupling parameter F*, at a fixed wave number k = 1.0. 

4-3. Quasi-two-dimensional Yukawa liquids confined by a parabolic potential 

The model adopted for the 2D Yukawa system, which assumes that the particles 
are constrained to move entirely within an ideal plane can be extended to describe 
more accurately the situation found in physical systems, by allowing small amplitude 
displacements of the particles perpendicular to the plane. In this extended model one 
applies a parabolic potential along the direction perpendicular to the plane, which 
then results in a quasi-two-dimensional confinement. Such confinement gives rise to a 
particle layer with finite width, or - at weaker confinement - to a sequence of multilayer 
structures, when the confinement or interaction potential is varied. The structural 
phase transitions (a change in the number of layers and in the accompanying crystal 
structures), relevant to particle traps, have theoretically been studied by Dubin [105j . 
while Totsuji et al. [106j . Bystrenko [107] as well as Qiao and Hyde |108j investigated the 
formation of layers in Coulomb and Yukawa systems in confined quasi-2D configurations. 

The number of layers formed in the liquid phase depends on the strength of the 
confinement. In contrast to the idealized 2D systems the layers have a finite width. Here 
we deal with the domain of parameters when a single layer is formed. At higher number 
of layers the mode structure is expected to be more complex |109l IllOl llllj . but the 
study of these modes is not within the scope of the present analysis. In a single-layer 
configuration the third degree of freedom of the particles, the displacement perpendicular 
to the plane, gives rise to an additional collective excitation, the "out-of-plane" V mode, 
besides the "in-plane" C and T modes found in (ideal) 2D layers. The out-of-plane mode 
in the crystallized state has been studied through simulations by Qiao and Hyde |112j . 
Results pertaining the strongly coupled liquid phase were analyzed in |113j . and will be 
summarized below. It should also be noted that a somewhat similar physical situation 
arises when a ID chain of particles is confined in the transverse direction by a parabolic 
potential: indeed, there is a similarity between the modes that represent excursions 
along the direction of the confining force in the ID and 2D systems. 

In the quasi-2D liquid system the particles can freely move in the (x, y) plane while 
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a confinement potential Vc{z) oc z"^ acts upon them when they are displaced from the z 
= plane. The confinement force is linear with respect to the "vertical" displacement, 

where the strength /q (besides F and k) is the third characteristic parameter of the 
system. At /o = 1 the confinement force at a vertical displacement 2; = a is equal 
to the magnitude of the force between two particles separated by a [defined by ([5])], 
interacting via Coulomb potential. Information about the (thermally excited) collective 
modes and their dispersion is obtained from the analysis of the correlation spectra of 
the longitudinal and (in-plane as well as out-of-plane) transverse current fluctuations. 
For the "in-plane" C and T modes we use eq. ( ITTi) . while for the out-of-plane mode the 
corresponding microscopic current 7r(fc,t) (which characterizes the V mode) is obtained 
as: 

7i{k,t) = k Vjzjt) exp[ikxj(t)]. (59) 
j 

Representative current fluctuation spectra for the three (£, V, and T) modes are 
displayed in figure [261 fo^' T = 100, /o = 2.0 and two different values of the screening 
parameter k = 0.27 and k, = 1.33. The frequency is normalized according to ([7]). We 
observe sharp peaks in the L{k,uj) spectra, similarly to the case of (ideal) 2D Coulomb 
and Yukawa liquids [521 [71], characteristic of long-lifetime collective excitations [114j . 
Peaks in the T mode spectra [see figure [26]^g,h)] show up only above a certain (cutoff) 
wave number, similarly to the case of 2D and 3D Yukawa systems, as discussed in the 
previous sections. 

The V mode possesses a finite frequency at = 0, which is, in general, characteristic 
of an optical mode. The first identification of this pseudo-optical behavior in a confined 
2D system is due to [114 ]. At small wave numbers the peaks of the spectra shift to lower 
a; as A; is increased. The width of the peaks of the P{k,uj) spectra become gradually 
broader when k is increased, as it can be seen in figure [261(e) and (f). It is noted that, on 
the other hand, the peaks become narrower as the strength of the confining potential, /o, 
is increased, which is an indication of an increasing lifetime of this collective excitation. 

The dispersion relations derived from the spectra are displayed in figure [27] for 
different values of /o and k for F = 100. At constant k, as shown in figure [27[(a). 
the frequency of the out-of-plane mode changes significantly as the strength of the 
confinement force, /o, is varied. The C and T modes are only slightly affected by the 
value of /q. The frequency of these modes is somewhat smaller at /q = 1.4, which is near 
the lower bound of /o for the formation of a single layer |113j . It is noted that at lower 
/o values, when two layers are formed, two longitudinal and two in-plane transverse 
modes appear, similarly to those identified in the classical (ideal) bilayer system |109j . 
Additionally, two out-of-plane transverse modes also emerge in the two-layered system, 
which are also believed to be in-phase and out-of-phase modes (when particles in the 
two layers oscillate in phase or with a phase difference of 180° in the two layers). The 
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Figure 26. Quasi-2D Yukawa liquids: dynamical structure function [S{k,uj)], 
longitudinal [L(fc,w)] and out-of-plane as well as in-plane transverse [P{k,uj) and 
T{k, uj)] current fluctuation spectra for k = 0.27 (a,c,e,g) and k= 1.33 (b,d,f,li). The S, 
L, and P spectra are plotted for the multiples of the smallest accessible wave number 
fcmin — 0.0886, while the T spectra are shown for higher wave numbers indicated by 
the labels in (g) and (h). The arrows in (a)-(f) indicate increasing wave numbers. 
Confining force: /q — 2. Reproduced from Ref. |113j . Copyright (2004) by the 
American Physical Society. 



Dynamical correlations and collective excitations of Yukawa liquids 



45 




Figure 27. Quasi-2D Yukawa liquids: dispersion relations for (a) k — 0.27 and 
different values of the amplitude /o of the confining potential, and (b) for fixed /g 
= 2 and different values of the screening parameter. Reproduced from Ref. [113j . 
Copyright (2004) by the American Physical Society. 
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Figure 28. Quasi-2D Yukawa liquids: the relation between the frequency oj{k — 0) 
and the Einstein frequency of the V mode for different values of the screening parameter 
K. Reproduced from Ref. |f 13j . Copyright (2004) by the American Physical Society. 

C mode exhibits a quasi acoustic behavior, with a hnear portion of the dispersion curve 
around k = 0, which widens with increasing k, as it can be seen in figure [271(b). The 
T mode shows an acoustic, ~ A; dispersion at small k, with a cutoff at a finite wave 
number. 

For the V mode du/dk < in the k < 2.1 domain. At higher wave numbers, the 
frequency of the mode slightly increases with k. This observation on the liquid system 
agrees well with that on the crystallized system |112] . where the same behavior was 
found, except that in the latter system the critical wave number (at which the group 
velocity duj/dk changes from negative to positive) also depends on the direction of the 
propagation. 

At A; = the whole layer oscillates in unison in the potential well with a frequency: 
= V%r2- (60) 
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A smaller confinement force results in a smaller uj{k = 0) and uj{k oo). At a constant 
/o the value of u{k = 0) does not change when k is varied, but - as shown in figure [271(b) 
- u!{k > 0) increases with decreasing k. This is explained by the decreased interparticle 
force (at an average particle separation) at higher k. 

The frequency of the out-of-plane mode at the A; — > oo limit (i.e. the Einstein 
frequency [351 ES]) can be calculated by considering the forces acting upon a single 
particle displaced in the z-direction, while all other particles are in rest in the z = 
plane. The force is the sum of the confining force and the force due to repulsion by the 
other particles, 

F{z) = -fo-^^z + F,{z). (61) 

The Fj.{z) contribution can be calculated as F^{z) = —dV^/dz, where Vj-^z) is the 
potential distribution due to a charge distribution p(x, y) in the z = plane. To obtain 
p{x, y) one may either use the radial (2D) pair correlation function (PCF) or consider 
the particles occupying hexagonal lattice sites in the z = plane. The Ft.{z) force is 
found to be a nearly linear function of the displacement z, in the \z\ < 0.3a domain, 
where the particle displacement is expected to fall. The resulting (Einstein) frequency 
(when the particles in the z = plane are situated at lattice sites) is [113] : 

^ ^{k oo) ^ / /o - 1.63exp(-1.37K) 
too ~ ioo ~\ 2 ■ ^ ' 

In the case of using the disordered configuration in the z = plane instead of lattice sites 
(through PCFs obtained in the liquid state simulations), a frequency very close to that 
given by fl62l) is obtained. At low values of k the Einstein frequency uje is significantly 
lower than uj{k = 0), as illustrated in figure l28l In the high k, limit the two frequencies 
are equal, as the screening becomes very strong and the particles interact very weakly. 
In this case the frequency of the V mode becomes nearly independent of k. 



5. Experimental results 

Experimental results on wave propagation and collective excitations in Yukawa systems 
have been accumulating in complex (dusty) plasma experiments since the mid 1990-s. 
Experiments have been carried out both on spontaneously generated and on externally 
excited waves. An early laboratory observation of longitudinal modes was reported by 
Barkan et al. |115j in 1995, followed by a more detailed study in 1997 |116j : these authors 
observed spontaneously generated waves in dust that filled a volume with a cyhndrical 
geometry. These waves grew as the result of the dust-acoustic instability, which was 
driven by an ion flow, and the experimenters were able to measure the wavelength 
and propagation speed. An early effort to excite waves by manipulation using an 
electrically-biased wire was reported by Pieper and Goree [117] . who also introduced 
a method of data analysis that yields the real and imaginary parts of the wave number, 
for the applied frequency. Repeating the measurements at various frequencies yielded a 
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dispersion relation. In some of these early experiments the ambient pressure was kept 
high in order to avoid instabilities. As a result, as pointed out by Rosenberg and Kalman 
[72], the waves were strongly damped, primarily by grain-neutral collisions. Thus, 
even though the experiments were conducted under strongly coupled conditions in the 
liquid state, the strong collisional damping washed away the difference between weakly 
coupled and strongly coupled dispersions (see [611 [72] for a more detailed discussion). 
Two recent experiments have further corroborated this picture. Bandyopadhyay et 
al. [11 8J investigated the acoustic dispersion of the longitudinal mode in the strong 
coupling regime over a wide range of the neutral pressure values and found du/dk < 
behavior of the dispersion curve that in the low collisional domain could be attributed 
to correlational effects. On the other hand, the experiment reported by Annibaldi et al. 
[119] . in the high collisional regime confirmed that in this domain the strong coupling 
effects were washed away completely. 

The first experiments where strong coupling effects were clearly displayed were 
done on a ID complex plasma in the crystalline state, realized as a chain of grains 
held together in the transverse direction by a confining potential. Longitudinal waves 
(along the direction of the chain) excited by the radiation pressure of a laser beam 
[120[ I121j in a parallel plate radio frequency discharge were observed: the analysis led 
to the conclusion that the weakly coupled theory of the longitudinal waves (referred 
to as "dust acoustic waves" [122]) was inadequate, while the description in terms of 
harmonic phonons of a system with short range interaction (referred to as "dust lattice 
waves" [123] ) provided a more satisfactory agreement with experiments. 

Generation of 3D complex plasmas in the laboratory under strong coupling 
conditions and at sufficiently low pressure, so that strong coupling effects become 
manifest turned out to be difficult. A good summary of the state of affairs as of 2000 is 
given by [64| . This paper and a later work |124j also discuss the spontaneous excitation 
of shear wave-like structures in a strongly coupled liquid at a low pressure. The plasma 
originally was in a 3D configuration, but assumed a layer structure in the course of the 
experiment, with particle excursions in the direction perpendicular to the layers. Thus 
it seems difficult to judge whether the observed waves were indeed shear waves or some 
more intricate excitation in the coupled layer system. 

The presence of the ion beam traversing the dust plasma generated in low- 
pressure gas discharges can create issues that can not be approached within the Yukawa 
model. Because of the anisotropy introduced by the ion beam, the Yukawa interaction 
will be modified in the vertical direction (along the beam), but probably not too 
much in the horizontal plane [125[ 1126] . This scenario was supported by experiment 
[127[ I128[ 11291 1130] . A further major problem due to the ion beam in the 3D geometry 
was identified by Joyce et al. |131] . that in the low pressure domain, where collective 
modes could be observable, ion-dust instability may lead to melting. To avoid these 
problems, most of the subsequent laboratory experiments were to favor 2D geometries 
over three-dimensional ones: in a 2D system these problems should be absent. Since 
1998 substantial progress in the understanding of the excitation and propagation of 
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waves in 2D Yukawa systems has ensued. In the strong couphng regime the constituent 
grains are, in principle, either in the crystaUine sohd or in the hquid state. In fact, 
in addition to the formation of large scale ordered lattice structures a more common 
configuration is an aggregate of micro-crystals whose prevailing disorder is expected to 
make the behavior of the aggregate quite similar to that of the liquid state. 

Longitudinal waves in a 2D dust plasma crystal were first observed experimentally 
in a parallel plate radio frequency discharge by Homann et al. |132] . The observation 
of transverse (in-plane) shear waves, the hallmark of strong coupling, excited by a 
chopped laser beam was reported by Nunomura et al. [85]. Their measurements of the 
dispersion relation revealed an acoustic, i.e., non- dispersive, character over the entire 
range of wave numbers measured, (0.3 < k < 1.2), at k ^ 0.74, with transverse sound 
speed and Einstein frequency values in agreement with theory |79l [80] . 

A series of beautiful experiments on the generation of Mach cones in the wake 
of an object moving through a 2D dusty plasma crystal was also crucial albeit in an 
indirect way, in determining the strong coupling characteristics of these systems. It 
was unambiguously shown |133l 11341 1135] that Mach cones appear when the velocity of 
the moving object (particle or laser spot) exceeds the longitudinal sound speed sl in 
the medium. Subsequent observations [136] with object velocities below this limit, but 
above the transverse shear sound speed st (the ratio of the two speeds had the value 
sl/st = 4.48 in the experiment) also demonstrated the excitation of a small angle Mach 
cone sustained by the transverse mode. 

More recent experiments done both on the 2D solid and liquid phases have been able 
to determine plasma parameters with sufficient accuracy and to perform measurements 
of great number of observables, so that detailed quantitative comparisons with the 
theoretical predictions have become possible. Experiments by Nunomura et al. |137] 
introduced laser manipulation, which avoided technical problems caused earlier by the 
electrical wires, and the dispersion relations were measured with greater accuracy. 
Subsequently Nunomura et al. [138] detected the spectra of self-generated longitudinal 
and transverse excitations along the two principal axes of a triangular lattice. The energy 
was concentrated along a well defined a;(/c) curve, representing the measured dispersion 
relation. The data covered the < k < 3.3 domain with k ^ 0.74; our comparison with 
the theoretical dispersion curves calculated by Peeters and Wu [HQ] and by Sullivan et 
al. [83] shows excellent agreement with these data. A partial frequency spectrum (i.e. 
density of states), based on the kinetic energy contents of the 4 selected modes was also 
generated: while comparison with the calculated spectrum (see section [3]) is possible, 
agreement beyond what is visible in figure [29] is not expected, since the theoretical 
spectrum includes all propagating modes [139] . In a subsequent work [140] the spectrum 
of waves was measured also for a number of directions in between the principal axes 
and over a much broader domain of wave number values: < < 6.6. In addition to 
the frequencies, the polarization angles of the modes were also determined (the mode 
polarizations can be described as "longitudinal" and "transverse" for propagation along 
the principal directions only). All these data show excellent agreement with theory 
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Figure 29. 2D system: density of states as obtained from the experiment of Nmiomura 
et al. jl38| and by theory. 




Figure 30. (color onhne) 2D system: wave dispersion in directions (a) 0, (b) 10, (c) 
20, and (d) 30 degrees (measured from the nearest neighbor direction) as obtained in 
the experiment of Zhdanov et al. |140j and calculated from lattice summation (heavy 
lines). Experimental data reproduced from Ref. |140j with kind permission of the 
authors. Copyright (2003) by the American Physical Society. 



Following a different line of approach Melzer )141] studied the normal modes of 
small 2D clusters of grains; what is of interest in the present context is the transition 
from the mode spectrum of a finite number of particles into that of an "infinite" 
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system. With a somewhat arbitrary assignment of labels for the normal modes, it 
was found that the average frequency as a function of k provides a fair resemblance to 
the Lj{k) dispersion of the infinite lattice, already for a cluster as small as consisting of 
34 particles. The scatter of frequencies around the uj{k) curve is of course substantially 
higher than in the experiments quoted above. (For somewhat related results see |34j). 
The theoretical understanding of the distribution of dynamical frequencies (as a function 
of temperature and particle number) presents a major theoretical challenge, with a 
very limited body of antecedents available in the literature [66l I142j . The possibility of 
generating experimentally observable scenarios from which information on the frequency 
distribution can be extracted should provide stimulus for new theoretical efforts. 

As to the liquid state, observational data available at the present time are quite 
recent and still rather limited. Nunomura et al. [15] studied the change of the thermally 
excited mode structure as the crystal lattice was melted and the system transited to the 
strongly coupled liquid state. The melting was achieved by directed laser heating. The 
theoretically predicted trends, such as the development of a cut-off wave number for the 
shear mode and the shift of the longitudinal mode frequency towards higher values are 
in fair agreement with MD data. On the other hand, the widths of the spectra seem 
to be higher than expected. It is difficult to relate the results of the MD studies of the 
break-up and isotropization of the lattice modes [82] to this experiment, since the F 
values where the observations of the liquid state were done are quite far from the phase 
transition point. 

A careful study of the transverse modes in the strongly coupled liquid state, in 
the vicinity of the melting point is due to Piel et al. [16j. These authors analyzed the 
propagation of externally excited shear waves, through a sophisticated data analysis 
technique that made it possible to collect information from a high noise environment. 
The experimental situation corresponded to k ~ 0.4, which allowed the comparison 
with the QLCA data for k = and k ^ 1 as lower and upper bounds. The authors 
found that within the wave number domain investigated (0 < A; < 2.5) the overall 
agreement between experiment and the QLCA model is quite satisfying. They note that 
even in the solid state the waves assume characteristics resembling those in the liquid 
state (angularly averaged dispersion), because the plasma crystal consists of domains of 
different orientations. For this reason there does not seem to be too much change in the 
dispersion, as one passes from the solid phase to the liquid phase. On the other hand, 
the damping is substantially higher on the liquid side and becomes stronger for low k 
values. 

Since the homogeneous liquid cannot sustain shear, the shear mode must vanish 
below some finite k^ value. The value of this cut-off wave number was recently studied 
by Nosenko et al. [17] in a low pressure experiment. At the k ~ 0.43 of the experiment 
kc ranges between 0.16 and 0.31, as the coupling strength F is varied from the melting 
value F = 155 down to F = 60. These values compare favorably with the values obtained 
by the MD simulations reported in [71] (see fig. [3Til . This can be taken as an indication 
that the cut-off is attributable to the intrinsic dynamics of the grains. Thus one can 
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conclude that at low pressures the contribution of the grain-neutral collisions to the 
generation of the cut-off is quite negligible. 




Figure 31. 2D Yukawa liquids: Transverse mode cutoff wavenumber fee. Experimental 
values are taken from [17 and are compared with 2D Yukawa molecular dynamics 
results. 



There are experiments on ID chains of grains held together in the transverse 
direction by a confining potential which reveal (in addition to the observation of 
longitudinal waves on such systems quoted above [12011121] ) excitations in the direction 
perpendicular to the axis of the chain. These modes bear physical features similar 
to the shear-waves in 2D systems. In particular, they exhibit the benchmark pseudo- 
optic behavior and the ensuing negative dispersion predicted by theory for the latter 
|108[ 11121 1113] . Experimentally, the mode dispersion was determined by analyzing 
spontaneously excited waves by Misawa et al. |143j and by creating the transverse 
waves through the manipulation of a single particle by Liu et al. |144l 1145] confirming 
these predicted features. In the 2D geometry, self-excited out-of-plane oscillations of 
particles were identified first by Nunomura [ 146] . Samsonov et al. |147] investigated 
the propagation of wave packets in the vertical (i.e. along the confinement) direction 
(see Section for more details) and confirmed the predicted dispersion characteristics 
[1081 1112[ 1113] of the out-of-plane (P, pseudo-optic) mode. 



6. Summary 

The objective of this review has been to summarize the huge body of information that 
has been gathered since 1990-s through theoretical analysis, computer simulations and 
laboratory experiments on the collective excitations of dusty (complex) plasmas and 
from this to determine the collective behavior of two- and three-dimensional strongly 
coupled Yukawa systems. The Yukawa model allows the mathematical analysis of an 
idealized system that represents a variety of actual many-particle physical systems 
(dusty plasmas, charged colloids, mesoscopic particles, etc.), which are characterized 
by (i) a significant ratio of the potential energy (originating from the high charge value 
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of the particles) to the kinetic energy in the system, as expressed through the plasma 
coupling parameter F, and (ii) by a particle-particle interaction that is strongly affected 
by a polarizable background coexisting with the main plasma. 

Two techniques have been used for the mathematical modeling: molecular 
dynamics, as a computational simulation method and the Quasi-Localized Charge 
Approximation as a theoretical scheme. The results generated by the two independent 
approaches have been found to be in excellent agreement with each other, and have 
been convincingly supported by the findings of laboratory experiments. Thus all these 
assembled data converge into a coherent and fairly complete physical picture, which 
has been presented in this Review. Nevertheless, there are quite a few areas that the 
reader may have expected to see in this paper, but which have been excluded from 
consideration. Thus some qualifying comments along this line are in order. 

• We have considered infinitely large, unbounded 2- and 3-dimensional systems only: 
thus effects relating to ID geometry, boundary conditions, inhomogeneities have 
been excluded. Some of these (Yukawa balls and disks e.g. [1481 1149"] ) have been 
attracting much attention lately. 

• A special configuration, familiar in semiconductor physics, is the bilayer geometry 
(consisting of two parallel 2D planes, separated from each other by a small distance 
d). While semiconductor devices are governed by Coulomb interaction, a similar 
configuration is of interest with systems where Yukawa interaction prevails [150j . 
The likelihood of the realization of such a geometry in laboratory dusty plasma 
experiments is not promising using identical grains, but may be more feasible in 
a microgravity environment. However, combining two species of differently sized 
microparticles in a conventional laboratory sheath geometry setup leads to the 
automatic formation of a bilayer configuration, due to the different Z/ m ratios, as 
recently pointed out by Matthews et al. |151j and demonstrated in a subsequent 
experiment by Smith et al. |152j . Bilayer systems posses rich variety of structural 
phases |150l I153[ I154j and a rather complex collective mode structure whose details 
exhibit a remarkable sensitivity to the layer separation |109[ I15H I155j . Further 
experimental investigation of this behavior would be of great interest. 

• Only single component systems (the Yukawa equivalent of the OCP, one 
component plasma) have been treated; the crucially important extension to two- or 
multicomponent cases (different masses, different charges) is not here. Creation of 
such systems in the current laboratory set-ups is hampered by technical reasons (but 
again may become feasible in a microgravity environment) and serious theoretical 
studies are lacking. 

• While we have presented simulation data of the dynamical fluctuation spectra in 
great detail, the evaluation and theoretical analysis of most of these data is still 
to be carried out. We have focused on the positions of the peaks in the spectra: 
the most important question amongst those whose analysis is incomplete is that of 
the widths of the peaks, which, in turn, are characteristic of the damping of the 
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excitations. What is missing primarily is a solid theoretical foundation through 
which the different mechanisms that lead to the damping of the collective modes 
could be reconciled and against which the simulation data could be tested. The 
QLCA analysis points at the main physical effects where the source of the damping 
should be sought, but no reliable analytic tool has emerged that would predict how 
the damping depends on the coupling strength and on the wavelength of the mode. 

• The theoretical tool (the QLCA) described in this Review is geared to strong 
coupling and it provides no linkage from the localization dominated strongly 
correlated behavior to the fluid-like weakly correlated behavior of the collective 
modes. Only more simulation work and a different theoretical approach would 
bridge this gap. 

• There are both some experimental |140j and simulation (see figure [8] and [32]) results 
available on the effect of the disorder on the collective mode spectrum. Both these 
and theoretical considerations suggest that one should think in terms of frequency 
distributions, rather than in terms of well-defined collective mode dispersions. 

• We have not discussed effects and phenomena relating to external or internally 
generated magnetic fields. These issues may become the topics of investigation 
for the next generation of complex plasma experiments. An externally imposed 
magnetic field could affect the polarizable medium (electrons and ions) and thus 
alter the effective interaction potential; at sufficient strength it may even change the 
orbits of the mesoscopic plasma particles and thus the prevailing mode structure 
|156[ 1157] . (As an example, consider a plasma with grains of i? = 1 micron, 
mass density p ~ 1.5 g/cm^ and Z ~ 3000; here a magnetic field of 2 Teslas 
would produce a dust cyclotron frequency Ucd ~ 0.16 rad/s, and with Td ~ 0.03 
eV a gyroradius r^yj-o ~ 0.5 cm. Thus, based on the criterion rgyro < confinement 
length, the creation of a magnetized plasma may become feasible. A more restrictive 
criterion may, however, emerge from the requirement Ucd > ^'coih the grain-neutral 
collision frequency. In a different vein, systems containing magnetically polarizable 
plasma particles are expected to exhibit a series of novel physical phenomena, both 
in equilibrium [158] and in terms of collective excitations [159] . 

• Transport coefficients may have been a legitimate topic for consideration in this 
Review, but partly for reason of economy, partly because their treatment requires 
a different (theoretical and simulation) methodology from those appropriate for 
the study of wave phenomena, the subject has not been included. The transport 
coefficients of 3D Yukawa systems in the liquid phase are relative well known. The 
self-diffusion was studied in [Mj, estimates for the viscosity were given in [160] . 
Molecular dynamics simulations have proven to be invaluable tools for studies of 
transport phenomena and made possible the determination of shear viscosity and 
thermal conductivity [HH 11611 11621 11631 11641 1165] . Recent theoretical work on 
this topic has focused on the mapping between Yukawa, Coulomb and hard-sphere 
systems |166[ 1167] . The effect of Langevin dynamics on the viscosity of 3D Yukawa 
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systems has been studied in |168j . For recent experimental work on 3D systems see 
e.g. [169]. 

The reahzation of 2D complex plasma liquids and the development of modern 
experimental (perturbation and data acquisition) techniques resulted in renewed 
interest of transport properties (which are especially interesting due to the 
controversies about the very existence of transport coefficients in low-dimensional 
systems). During the last few years several experimental and simulation studies 
have been carried out on the shear viscosity [131 El 11701 11711 I172j . thermal 
conductivity |173l 1174] and diffusion [100[ llOlj properties of 2D Yukawa liquids, 
and this topic is expected to attract further attention. 
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